Bike Sharing

Introduction

• Keep in mind the tradeoffs based on speed, explainability, simplicity & performance
• Bike sharing data is found on Kaggle
• Potential things to explore:
  • Duration of travel, departure and arrival position
  • Casual vs. registered rental wrt weather, time, etc.
  • Understand seasonality / time series analysis
  • Predict number of rentals on a day
  • Event and Anomaly Detection

---

Data

Sample

options(scipen=999)
seed = 2018
pkgs = c('data.table', 'magrittr', 'caret', 'glmnet', 'ranger', 'plotly', 'iml')
inst = lapply(pkgs, library, character.only=T, quietly=T); rm(pkgs, inst)
source('~/Documents/R/Utils/functions/summarize.R')
source('~/Documents/R/Utils/functions/one_hot.R')
source('~/Documents/R/Utils/functions/model_auc.R')
source('~/Documents/R/Utils/functions/feature_comparison.R')
source('~/Documents/R/Utils/functions/feature_distribution.R')
source('~/Documents/R/Utils/functions/helper.R')
day = fread("~/Rrrrr/projects/bike_sharing/day.csv") %>%
  .[, dteday := as.Date(dteday)]
hour = fread("~/Rrrrr/projects/bike_sharing/hour.csv") %>%
  .[, dteday := as.Date(dteday)]
head(day)

---

Summarize

The goal is to find variables with:

• missing values
• low variance
• outliers – casual users shows the most diverse behavior

Day

sum.day = Summarize(day); DT::datatable(sum.day, options = list(dom = "tip"))

---

Hour

sum.hour = Summarize(hour); DT::datatable(sum.hour, options = list(dom = "tip"))

---

Learnings

• Holiday is only marked for weekdays

day[,.N,.(weekday, workingday, holiday)]

---

Anomaly detection EDA

• Given condition, actual rental is way higher than predicted rentals, especially casual renters on days of events
• Time of rental – might be a peak at certain time
• My guess is holiday is more likely to have events – only 21 holidays
• Holiday & weekends people might use bikes on a longer stretch of time
• Clustering on raw features & on shapley values

Feature Density

FeatureComparison(list(holiday = day[holiday==1, cnt], 
                       weekday = day[holiday==0 & workingday==1, cnt],
                       weekend = day[holiday==0 & workingday==0, cnt]),
                  vectorName = "Total rental", sample = F, numerical = T)

---

FeatureComparison(list(holiday = day[holiday==1, casual], 
                       weekday = day[holiday==0 & workingday==1, casual],
                       weekend = day[holiday==0 & workingday==0, casual]),
                  vectorName = "Casual rental", sample = F, numerical = T)

---

FeatureComparison(list(holiday = day[holiday==1, registered], 
                       weekday = day[holiday==0 & workingday==1, registered],
                       weekend = day[holiday==0 & workingday==0, registered]),
                  vectorName = "Registered rental", sample = T, numerical = T)

Feature Over Time

• 2012-10-29 & 2012-10-30 why is it so low?
  • Hurricane Sandy

FeatureDistribution(list(weekday = day[holiday==0 & workingday==1],
                         weekend = day[holiday==0 & workingday==0],
                         holiday = day[holiday==1], overall = day),
                    featureName = "cnt", periodName = "dteday")

---

• Casual users show the most diverse behaviors between holiday, weekends & weekdays
• 2011-07-04 is very high – 4th of July!

FeatureDistribution(list(weekday = day[holiday==0 & workingday==1],
                         weekend = day[holiday==0 & workingday==0],
                         holiday = day[holiday==1], overall = day),
                    featureName = "casual", periodName = "dteday")

---

• Looks like registered users like to use more during weekdays, which makes sense as the registered users use bikes as work commute transportation

FeatureDistribution(list(weekday = day[holiday==0 & workingday==1],
                         weekend = day[holiday==0 & workingday==0],
                         holiday = day[holiday==1], overall = day),
                    featureName = "registered", periodName = "dteday")

---

Weather

• Casual users like to bike on cooler days during summer, and clearer days in spring & fall

FeatureDistribution(list(Clear = day[weathersit==1], 
                         Mist = day[weathersit==2],
                         Light = day[weathersit==3]),
                    featureName = "casual", periodName = "mnth")

---

• Again, weather has less effect on registered users (similar to weekdays)

FeatureDistribution(list(Clear = day[weathersit==1], 
                         Mist = day[weathersit==2],
                         Light = day[weathersit==3]),
                    featureName = "registered", periodName = "mnth")

---

Day of week

• Casual rental: Saturdays & Sundays are heavy days

FeatureComparison(list(S = day[weekday==0, casual], 
                       M = day[weekday==1, casual],
                       T = day[weekday==2, casual],
                       W = day[weekday==3, casual],
                       Th = day[weekday==4, casual],
                       F = day[weekday==5, casual],
                       St = day[weekday==6, casual]),
                  vectorName = "Casual rental", sample = F, numerical = T)

---

• Registered rental: Wednesdays & Thursdays are heavy days

FeatureComparison(list(S = day[weekday==0, registered], 
                       M = day[weekday==1, registered],
                       T = day[weekday==2, registered],
                       W = day[weekday==3, registered],
                       Th = day[weekday==4, registered],
                       F = day[weekday==5, registered],
                       St = day[weekday==6, registered]),
                  vectorName = "Registered rental", sample = F, numerical = T)

---

Time

• Starting from 0, every 15\(^\circ\) is one hour; so 120\(^\circ\) means 8am and 255\(^\circ\) means 5pm

• How far the dot is from origin is the average number of users during that time

• Casual users
  • Rentals on Saturday > Sunday > Weekdays
  • More rentals in the afternoon on weekends
  • Most weekday rental at ~5pm
• Registered users
  • More rentals during weekdays > weekends, Saturday & Sunday uses are not very different
  • Peak hours are 8am, 5pm & 6pm during weekdays

time = hour[workingday==1, .(day = "Weekday", c = mean(casual), r = mean(registered)), hr] %>%
  rbind(hour[weekday==6, .(day = "Sat", c = mean(casual), r = mean(registered)), hr]) %>%
  rbind(hour[weekday==0, .(day = "Sun", c = mean(casual), r = mean(registered)), hr]) %>%
  .[, deg := hr*15] %>%
  .[, let := LETTERS[hr+1]]
# plot_ly(time, x=~hr, y=~cnt, color = ~day) 
p <- plot_ly(type = 'scatterpolar', mode = 'markers', 
             r = time[day=="Weekday"]$c,
             theta = time[day=="Weekday"]$deg,
             name = "Casual_Week", marker = list(opacity = 0.5)) %>%
  add_trace(type = 'scatterpolar', mode = 'markers', 
            r = time[day=="Sat"]$c,
            theta = time[day=="Sat"]$deg,
            name = "Casual_Sat", marker = list(opacity = 0.5)) %>%
  add_trace(type = 'scatterpolar', mode = 'markers', 
            r = time[day=="Sun"]$c,
            theta = time[day=="Sun"]$deg,
            name = "Casual_Sun", marker = list(opacity = 0.5)) %>%
  add_trace(type = 'scatterpolar', mode = 'markers', 
            r = time[day=="Weekday"]$r,
            theta = time[day=="Weekday"]$deg,
            name = "Reg_Week",
            marker = list(opacity = 0.5)) %>%
  add_trace(type = 'scatterpolar', mode = 'markers', 
            r = time[day=="Sat"]$r,
            theta = time[day=="Sat"]$deg,
            name = "Reg_Sat",
            marker = list(opacity = 0.5)) %>%
  add_trace(type = 'scatterpolar', mode = 'markers', 
            r = time[day=="Sun"]$r,
            theta = time[day=="Sun"]$deg,
            name = "Reg_Sun",
            marker = list(opacity = 0.5)) %>%
  layout(polar = list(radialaxis = list(angle = 0),
                      angularaxis = list(direction = 'clockwise', dtick = 15)))
p

---

Holiday

• 21 Holidays – easy to spot days with events

day[,.N, .(holiday, workingday)]

day[holiday==1,.(dteday, yr, weekday, casual, registered, cnt)]

---

• Time

hour[dteday=="2011-07-04",.(hr, weathersit, temp, atemp, hum, windspeed, casual, registered, cnt)]

---

Regression

Train-test split

• Time split

cutoff = "2012-09-01"
# train-test split by time
data.train = hour[dteday<cutoff, -c("instant", "dteday", "casual", "registered")]
week.train = hour[dteday<cutoff] %>% .[, week := week(dteday)] %>%
  .[, -c("instant", "dteday", "casual", "registered")]
date.train = hour[dteday<cutoff] %>% .[, week := week(dteday)] %>%
  .[, date := mday(dteday)] %>%
  .[, -c("instant", "dteday", "casual", "registered")]
cas.train = hour[dteday<cutoff]$casual
reg.train = hour[dteday<cutoff]$registered
data.test = hour[dteday>=cutoff, -c("instant", "dteday", "casual", "registered")]
week.test = hour[dteday>=cutoff] %>% .[, week := week(dteday)] %>%
  .[, -c("instant", "dteday", "casual", "registered")]
date.test = hour[dteday>=cutoff] %>% .[, week := week(dteday)] %>%
  .[, date := mday(dteday)] %>%
  .[, -c("instant", "dteday", "casual", "registered")]
# check to make sure distribution is consistent
FeatureDistribution(list(Train = hour[dteday<cutoff], 
                         Test = hour[dteday>=cutoff]),
                    featureName = "cnt", periodName = "dteday")

---

One-hot-encoding

• Season & weekday

hour.ohot = copy(hour)
OneHot(hour.ohot, "mnth")
OneHot(hour.ohot, "hr")
OneHot(hour.ohot, "season")
OneHot(hour.ohot, "weekday")
OneHot(hour.ohot, "weathersit")
ohot.train = hour.ohot[dteday<cutoff, -c("instant", "dteday", "casual", "registered")]
ohot.test = hour.ohot[dteday>=cutoff, -c("instant", "dteday", "casual", "registered")]

---

Additional work

• Identify & remove outliers

---

Model build

• Train with entire dataset
• Train with events (0/1) column
• Train separate models for registered vs. casual users
  • Exact
• Evaluation metric determines how good a model is, for e.g. RMSE means predicted values are in a narrow band in comparison to MAE

GLM

# 10-fold CV
set.seed(seed)
glm.ctr <- trainControl(method = "repeatedcv", number = 10, repeats = 5)
# search grid
glm.grd <-  expand.grid(alpha = 0:10/10, lambda = 10^seq(3,-6))
# model                        
set.seed(seed)
glm.all <- train(x = ohot.train[,-"cnt"], y = ohot.train$cnt, 
                 method = "glmnet", trControl = glm.ctr, tuneGrid = glm.grd)

There were missing values in resampled performance measures.

plot(glm.all)

---

glm.all$bestTune

---

• RMSE = 132.4112055
• RMSE results in conservative predictions – predicted values are all < 600

pred.glm = data.table(actual = ohot.test$cnt, 
                      pred = predict(glm.all, newdata = ohot.test))
plot_ly(pred.glm, x = ~actual, y = ~pred) %>%
  add_markers() %>%
  add_lines(x=c(1,800), y=c(1,800), showlegend = F)

---

• In comparison, RMSE on train = 95.1339412
• GLM has a higher bias

pred.glm.train = data.table(actual = ohot.train$cnt, 
                            pred = predict(glm.all, newdata = ohot.train))
plot_ly(pred.glm.train, x = ~actual, y = ~pred) %>%
  add_markers() %>%
  add_lines(x=c(1,800), y=c(1,800), showlegend = F)

---

GLM(2models)

Casual model

# model                        
set.seed(seed)
glm.cas <- train(x = ohot.train[,-"cnt"], y = cas.train, 
                 method = "glmnet", trControl = glm.ctr, tuneGrid = glm.grd)

There were missing values in resampled performance measures.

plot(glm.cas)

---

glm.cas$bestTune

---

Registered model

set.seed(seed)
glm.reg <- train(x = ohot.train[,-"cnt"], y = reg.train, 
                 method = "glmnet", trControl = glm.ctr, tuneGrid = glm.grd)

There were missing values in resampled performance measures.

plot(glm.reg)

---

glm.reg$bestTune

---

Performance

glm.two = data.table(hour[dteday>=cutoff, .(actual = cnt, dteday)],
                      cas = predict(glm.cas, newdata = data.fe)[14492:17379],
                      reg = predict(glm.reg, newdata = data.fe)[14492:17379]) %>%
  .[, pred := cas + reg]

Error in cbind2(1, newx) %*% nbeta : 
  Cholmod error 'X and/or Y have wrong dimensions' at file ../MatrixOps/cholmod_sdmult.c, line 90

---

FeatureComparison(list(glm.two$actual, glm.two$pred))

---

• Average hourly prediction

FeatureDistribution(list(Train = hour[dteday<cutoff,.(dteday, group = "train", value = cnt)],
                         Test = glm.two[,.(dteday, group = "actual", value = actual)],
                         Preds = glm.two[,.(dteday, group = "pred", value = pred)]), 
                    featureName = "value", periodName = "dteday")

---

XGB

• nrounds (# Boosting Iterations), max_depth (Max Tree Depth), eta (Shrinkage), gamma (Minimum Loss Reduction), colsample_bytree (Subsample Ratio of Columns), min_child_weight (Minimum Sum of Instance Weight), subsample

# 10-fold CV
set.seed(seed)
xgb.ctr <- trainControl(method = "repeatedcv", number = 10, repeats = 1)
# search grid
# glm.grd <-  expand.grid(alpha = 0:10/10, lambda = 10^seq(3,-6))
# model                        
set.seed(seed)
xgb.all <- train(x = data.train[,-"cnt"], y = data.train$cnt, 
                 method = "xgbTree", trControl = xgb.ctr)
plot(xgb.all)

---

xgb.all$bestTune

---

Performance

• RMSE = 79.3423604

pred.xgb = data.table(actual = data.test$cnt, 
                      pred = predict(xgb.all, newdata = data.test))
plot_ly(pred.xgb, x = ~actual, y = ~pred) %>%
  add_markers() %>%
  add_lines(x=c(1,800), y=c(1,800), showlegend = F)

---

• In comparison, RMSE on train = 45.8525654

pred.xgb.train = data.table(actual = data.train$cnt, 
                            pred = predict(xgb.all, newdata = data.train))
plot_ly(pred.xgb.train, x = ~actual, y = ~pred) %>%
  add_markers() %>%
  add_lines(x=c(1,800), y=c(1,800), showlegend = F)

---

XGB(2models)

• Slightly performance as the normal XGB model

Casual model

set.seed(seed)
xgb.cas <- train(x = data.train[,-"cnt"], y = cas.train, 
                 method = "xgbTree", trControl = xgb.ctr)
plot(xgb.cas)

---

xgb.cas$bestTune

---

Registered model

set.seed(seed)
xgb.reg <- train(x = data.train[,-"cnt"], y = reg.train, 
                 method = "xgbTree", trControl = xgb.ctr)
plot(xgb.reg)

---

xgb.reg$bestTune

---

Performance

• RMSE = 77.6327841

pred.two = data.table(hour[dteday>=cutoff, .(actual = cnt, dteday)],
                      cas = predict(xgb.cas, newdata = data.test),
                      reg = predict(xgb.reg, newdata = data.test)) %>%
  .[, pred := cas + reg]
plot_ly(pred.two, x = ~actual, y = ~pred) %>%
  add_markers() %>%
  add_lines(x=c(1,800), y=c(1,800), showlegend = F)

---

FeatureComparison(list(pred.two$actual, pred.two$pred))

---

• Average hourly prediction

FeatureDistribution(list(Train = hour[dteday<cutoff,.(dteday, group = "train", value = cnt)],
                         Test = pred.two[,.(dteday, group = "actual", value = actual)],
                         Preds = pred.two[,.(dteday, group = "pred", value = pred)]), 
                    featureName = "value", periodName = "dteday")

---

XGB(week)

# 10-fold CV
set.seed(seed)
xgb.ctr <- trainControl(method = "repeatedcv", number = 10, repeats = 1)
# search grid
# glm.grd <-  expand.grid(alpha = 0:10/10, lambda = 10^seq(3,-6))
# model                        
set.seed(seed)
xgb.week <- train(x = week.train[,-"cnt"], y = week.train$cnt, 
                 method = "xgbTree", trControl = xgb.ctr)
plot(xgb.week)

---

xgb.week$bestTune

---

• RMSE = 78.897122

week.xgb = data.table(hour[dteday>=cutoff, .(actual = cnt, dteday)],
                      pred = predict(xgb.week, newdata = week.test))
plot_ly(week.xgb, x = ~actual, y = ~pred) %>%
  add_markers() %>%
  add_lines(x=c(1,800), y=c(1,800), showlegend = F)

---

• Average hourly prediction

FeatureDistribution(list(Train = hour[dteday<cutoff,.(dteday, group = "train", value = cnt)],
                         Test = week.xgb[,.(dteday, group = "actual", value = actual)],
                         Preds = week.xgb[,.(dteday, group = "pred", value = pred)]), 
                    featureName = "value", periodName = "dteday")

---

XGB(date)

set.seed(seed)
xgb.date <- train(x = date.train[,-"cnt"], y = date.train$cnt, 
                 method = "xgbTree", trControl = xgb.ctr)
plot(xgb.date)

---

xgb.date$bestTune

---

• RMSE = 77.7874697

date.xgb = data.table(hour[dteday>=cutoff, .(actual = cnt, dteday)],
                      pred = predict(xgb.date, newdata = date.test))
plot_ly(date.xgb, x = ~actual, y = ~pred) %>%
  add_markers() %>%
  add_lines(x=c(1,800), y=c(1,800), showlegend = F)

---

• Average hourly prediction

FeatureDistribution(list(Train = hour[dteday<cutoff,.(dteday, group = "train", value = cnt)],
                         Test = date.xgb[,.(dteday, group = "actual", value = actual)],
                         Preds = date.xgb[,.(dteday, group = "pred", value = pred)]), 
                    featureName = "value", periodName = "dteday")

---

XGB(no mnth)

set.seed(seed)
xgb.mth <- train(x = data.train[,-c("mnth","cnt")], y = data.train$cnt, 
                 method = "xgbTree", trControl = xgb.ctr)
plot(xgb.mth)

---

xgb.mth$bestTune

---

Performance

• RMSE = 79.588598

mth.xgb = data.table(hour[dteday>=cutoff, .(actual = cnt, dteday)],
                     pred = predict(xgb.mth, newdata = data.test))
plot_ly(mth.xgb, x = ~actual, y = ~pred) %>%
  add_markers() %>%
  add_lines(x=c(1,800), y=c(1,800), showlegend = F)

---

• In comparison, RMSE on train = 46.5857025

mth.xgb.train = data.table(actual = data.train$cnt, 
                           pred = predict(xgb.mth, newdata = data.train))
plot_ly(mth.xgb.train, x = ~actual, y = ~pred) %>%
  add_markers() %>%
  add_lines(x=c(1,800), y=c(1,800), showlegend = F)

---

Performance

FeatureDistribution(list(Train = hour[dteday<cutoff,.(dteday, group = "train", value = cnt)],
                         Test = date.xgb[,.(dteday, group = "actual", value = actual)],
                         NoMth = mth.xgb[,.(dteday, group = "pred", value = pred)],
                         XGB2 = pred.two[,.(dteday, group = "pred", value = pred)],
                         Week = week.xgb[,.(dteday, group = "pred", value = pred)],
                         Date = date.xgb[,.(dteday, group = "pred", value = pred)]), 
                    featureName = "value", periodName = "dteday")

---

Explain

#calc.relimp

---

save.image(file = "bike.RData")
load("bike.RData")